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Separating Flow Behind a Back- Step. Part I. 
Quasi-Hydrodynamic Equations and Computation of a 

Laminar Flow 

We demonstrate the results of the numerical modelling of a plane two-dimensional viscous in- 
compressible flow in a channel with a back-step. As a mathematical model we take equations for a 
incompressible flow based on the quasi-hydrodynamic (QHD) equations. We present a phenomeno- 
logical derivation of the QHD equations and show a relation of these equations to the Navier-Stokes 
system. We test the proposed numerical algorithm by computing of certain laminar flows. 

1 Introduction 

The paper is devoted to the numerical modelling of the viscous incompressible flow behind the 
back-step in the channel with a sudden broadening. The quasi-hydrodynamic (QHD) system of 
equations is used as a mathematical model. £Q, 

In the second part of our paper we discuss the derivation of the QHD equations and their 
relation to the Navier-Stokes equations. The QHD equations broaden the possibilities of the 
classical Navier-Stokes model in description of the viscous compressible gas flows. When the 
Navier-Stokes equations are applicable, the additional dissipation of the QHD equations makes 
little influence on the solutions, but provides the stability of numerical computations. In certain 
cases of weakly rarified flows the QHD equations give the solution that describes experimental data 
better than the Navier-Stokes model 

Probation of the QHD equations for computing the incompressible liquid flows and for the 
problems of thermal and thermocapillar convection was carried out in @]-[7|. In particular, it was 
shown that the QHD equations are effective for modelling of nonstationary flows 6 . 

The size of the separation zone behind the step is a sensitive characteristic feature of laminar 
flows, which strongly depends on the flow velocity and on the geometry of the region considered 
in the problem. Analytical expression for the dependence of the separation zone's length on the 
Reynolds number and the relative height of the step for two-dimensional flows is given, for example, 
in 8 . The length of the separation zone grows almost linearly with the increasing of the Reynolds 
number. Laminar flows behind the back step are well simulated numerically and the results of 
two-dimensional computations of different authors are in a good agreement with experimental 
data [8]-|ll|. It allows us to use this problem as a test for probation of new numerical algorithms. 

In this paper we present the computations of the laminar flows in comparison with the previously 
published results. Data obtained from literature is used for evaluating the robustness and the 
accuracy of the numerical algorithm for computing the flow behind the step, which is based on the 
QHD equations. In the continuation of this paper (Part II) the method proposed here is applied 
to the numerical investigation of turbulent flows behind the step. 



2 Mathematical model 

In this section we describe the physical principles that form the basis for the phenomenological 
derivation of the new quasi-hydrodynamic (QHD)system of equations. Using time-space averaging 
for introducing the principal hydrodynamic values - density, velocity and temperature - is the 
essential and fundamental feature that distinguishes our method from the Navier-Stokes theory, 
where hydrodynamic values are introduced based on space averaging. 



2.1 Integral conservation laws 



Let us consider an inertial Cartesian coordinate system (x\, x 2 , x 3 ) in the Euclidian space i?| . 
Let (e*i, e*2, §3) be the corresponding orthonormal basis of unit vectors and let us denote time as t. 
We shall use the following standard notation for the variables describing the viscous compressible 
thermoconducting flow: p — p(x, t) - density, u = u(x, t) - velocity p = p(x, t) - pressure, 
£ = e(x, t) - specific internal energy, T = T(x,t) - temperature, s — s(x,t) - specific entropy. 

Suppose that the medium is two-parametric, that is, only two out of five thermodynamical 
parameters p, p, e, T, s are independent, and we are given by the state equations 

p=p(p,T), s = e(p,T), s = s(p,T). (1) 

Let F = F(x, t) be the mass density of external forces. For example, in case of the liquid in the 
gravitational field of the Earth it will be F — g, where g is the gravity acceleration. 
Our first postulate is the law of conservation of mass in the following form: 

dp 

— + divj m = 0. (2) 
at 

We suppose that the mass flux density vector j m — j m (x,t) is defined in every point x of the flow 
in every moment of time t. In the region occupied by the flow we take an arbitrary moving material 
volume V = V(t) with the smooth surface E = E(t), oriented with the field of external normal unit 
vectors n. We also suppose that the volume V(t) originates from the volume Vb = V(to), where to 
is the initial moment of time, by continuous deformation, caused by the motion of particles Vo along 
the trajectories, determined by the vector field j m /p- Using the well-known Euler-Liouville 
identity 

-J<pdV = J[D<p + <pdw(j m /p)]dV, (3) 

where tp — (p(x,t) is a certain continuously diffcrcntiable scalar or vector field, dV is a volume 
element in i?| and D — d/dt + {j m /p) • V is the differential operator, we present the law of 
conservation of mass (J2| in the integral form: 



d 



pdV = 0. (4) 



v 



The second postulate is the law of conservation of momentum 



j t J^dV = J pFdV + j "J{n- P)dE, (5) 

where c?E is the element of the surface E in the vicinity of the unit vector n. The rate of variation 
of the momentum in the volume V equals to the sum of all forces applied to it. The first integral 
in the right hand side of JSJ) is a volume force caused by the external field; the second stands 
for the forces, caused by pressure and internal viscous friction, that are applied to the surface E. 
The variable P = P(x,t) is called the tensor of internal tensions. The symbol (n ■ P) means the 
contraction (dot product) of the vector n and the second rank tensor P with respect to the first 
index of the tensor. Respectively, (P ■ n) means that the contraction of P and n is done with 
respect to the second index of P. If the tensor P is symmetric, then {ft - P) = (P ■ n). 
The third postulate is the law of conservation of the total energy 



Jt J p (T + £ ) dv = j <3™-FW + /£ (i: '" )dE ~/X ( ^" )dE - (6) 

Here the first integral in the right hand side of I© equals to the capacity of the external volume 
forces that are applied to the volume V; the second is understood as the capacity of the surface 
forces of the pressure and the internal viscous stress. The last term in I© describes the influx of 
energy in a single unit of time through the surface E due to the processes of the heat transfer. 
Actual expressions for the vector fields A — A(x, t) and q = q(x, t) will be given below. 



The fourth postulate expresses the law of conservation of the moment of momentum: 



[xx(pu)]dV = J [xxpF]dV + J J [x x (ft ■ P)]dE. (7) 

It is presented in its classical form. Internal moments and the distributed mass and surface pairs 
are not taken into consideration. The symbol x denotes the cross product of two vectors. 
Our fifth postulate is the second law of thermodynamics. It looks as follows: 

</-i</r = - / / / xdv. (8) 



dt ./,/ ./ .hp T 



V 



The surface integral in the right hand side JSJ defines the rate of the variation of entropy in the 
volume V due to the thermal flux. It may be both positive or negative. The last integral is always 
non-negative: it gives the production of entropy due to the internal irreversible processes. 



2.2 Transfer to differential equations 

Just like in the case of the Navier-Stokes system , for the transfer from the integral relations 
(0J-JHJ) to the corresponding differential ones we use the Liouville formula for differentiating 
the integral over the moving material volume. Doing it, we shall suppose that all the principal 
macroscopic parameters of the medium are sufficiently smooth functions of time and spatial coor- 
dinates. Taking in consideration that the volume V is arbitrary, we obtain differential equations 
for the balances of the mass 

dp 

— + divj m = 0, (9) 



of the momentum 



of the total energy 

d_ 

dt 



d(pu) 
dt 



dt 



+ div(j m (g> u) = pF + divP, 



(u 



div 



of the moment of momentum 



— [fx pu] + div{jr, 
at 



/ u 2 



— (jm ' F) + div A — divq, 



i) 



® [x x u]) = [x x pF] + -^[x x PijSji 



(10) 



(11) 



(12) 



and of the entropy 



d{ps) 
dt 



+ div(j m s) = -div{^j + X 



(13) 



Here ( jm ® u) is the second rank tensor obtained as a direct product of the vectors j m and it. When 
we take the divergence of the second rank tensor, we carry out the contraction with respect to its 
first index. The symbol in the equation i112[l means the portrait of the tensor P in the basis 
(e*i, e*2, ez). The summation is carried out with respect to the indexes i and j that appear twice. 

The system (f?Tll — (f 1 3|) is not closed. It is necessary to introduce the variables j m , P, q, A, X as 
the functions of macroscopic parameters of the medium and their derivatives. The closure problem 
can be solved in several ways. 



2.3 The classical approach to the closure problem. The Navier— Stokes 
equations 

First of all let us discuss the classical approach \T2\ , in which the averaging over a certain set of 
physically infinitely small volumes from the space i?| at the fixed moment of time t is used for 
definition of hydrodynamic variables. In this case the mass flow density vector j m at the arbitrary 



point (x, t) coincides with the average momentum of the unit volume pu, so the first closure relation 
looks as follows: 

jm = PU. (14) 

After that the pressure and inner viscous friction forces are introduced. They act instantly on the 
surface of the material volume. The law of motion of the latter is chosen in the same way as in 
the rigid body mechanics. In literature this assumption is called the solidification principle. The 
balance equation for the angular momentum i|12|) follows from the momentum conservation law 
<|1U|) and the symmetry of the tension tensor P. In the theory for Newtonian media P = Pns is 
defined by the expression 

p = n NS -pi, (is) 

where 

U NS = r)[(V ® u) + (V ® u) T - (2/3)Idivu\ (16) 

is the Navier-Stokes shear-stress tensor, I the unit tensor is the invariant of the second range. 
The heat flow q = qNS is defined according to the Fourier law 

q = -aiVT. (17) 

The hypothesis H16fl and l|17|) for the ideal monoatomic gases with small Knudsen numbers are 
confirmed by the kinetic computations. The work of the surface pressure forces and inner viscous 
shear-stress forces in a unit of time is computed using the same formula as in the rigid body 
mechanics, that is: 

A={P NS -u). (18) 
The specific thermodynamical entropy is supposed to satisfy the Gibbs differential identity 

Tds = de+pd(l/p). (19) 

Its balance equation l|13|) may be obtained as a consequence of the mass, momentum and energy 
conservation laws lflH|l - fTT|) . if we choose X = Xns as 

x=*(m 2 + inm: " Ns \ (20) 

V T J 2r)T ' v ; 

where (11^5 : Uns) — Y^i j=i(J^-Ns)ij(P-Ns)ij ~ is the double dot product of two identical tensors. 
Note that the right hand side of the equality (|20() is non-negative. Substitution of the expressions 
l|14 |) — (flS p into equations © - (|10|l gives us the classical Navier-Stokes system. The dependensies 
r] = 7y(p, T) and ae = ae(p, T) may be either found experimentally or derived from the kinetic theory 
of gases. 

2.4 The non-traditional approach to the closure problem. The quasi- 
hydrodynamic system 

Another way of solving the problem of closing the system (fTTft — (| IL 3 ft was proposed by Yu.V.Sheretov 
in 0, HI To define hydrodynamic variables he used not the spatial, but the time-spatial averaging 
over a certain set of physically infinitely small four-dimensional volumes in the space i?| t . He has 

proved that in the case of such time-spatial averaging the mass flow density vector j mi generally 
speaking, doesn't coincide with the average momentum of the unit volume pu. Detailed analysis 
of different possibilities for choosing the variables j m , P, q, A X gave the following result: 

j m = p(u~w), (21) 
P = -pi + Il N s + pu®w, (22) 



q = -aeVT, 



(23) 



A = (TLns • u) + pu{w ■ u) — p(u — w), 



(24) 



X = ac 



VT\2 (77 ws : n NS ) pw'< 



2r/T 



tT ' 



where 



w = — [p(u ■ V)u + Vp — pF] . 
P 



(25) 



(26) 



The parameter t = r(p, T) describes the scale of temporal smoothing. The formula for computing 
this parameter was proposed in [2] : 

= o 2' ^ 2? ) 

where 7 is the isentropic exponent, Sc is the Schmidt number (which is close to 1 for gases), c s is 
the sound velocity. The value of r agrees by order with the average mean free path of the particles 
in gas. Computations for moderately rarified gases confirm the correctness of this choice of the 
smoothing parameter 

Having substituted the expressions (|2T|) . and (PH instead of j m , P and A in i|§|l- ljTT|) . we 

obtain the quasi-hydrodynamic (QHD) system of equations: 



d(pu) 
dt 



— - + divlpu) = div(pw), 
at 



+ div{pu ® u) + X7p = pF + divllNs + div[(pw ® u) + (pu (g> w)}, 



(28) 
(29) 



d_ 



p^— + s^j + div P"u(^— + + pu + divq = pF ■ (u — w) + 

r / u \ 
div(II^s ■ u) + div pwy — + ej + pw + pu(w ■ u) 



(30) 



The QHD system I|28|) -(I30 |I becomes closed if it is equipped with the state equations (JIJ, and the 
coefficients 77, ae and r are presented as functions of macroscopic parameters of the media. The 
substitution of the expressions l|21|l , (|23[l and (|25|l into (|13H gives the entropy balance equation 



d(ps) 
dt 



div(pus) = div(pws) + div ^ae ^ + ae^ '— 



9. 



QHD 



(31) 



in which 



QHD 



(Hns ■ n^g) 

2^ 



pw 2 



is the non-negative dissipative function. 

A number of theoretical results is obtained for the QHD system <|28(1 - <|3U[) m CD- (2- I n par- 
ticular, it has been shown that the stationary QHD system in dimensionless variables differs from 
the corresponding Navier-Stokes equations only in terms of the second order of magnitude with 
respect to the Knudsen number. Its laminar boundary layer approximation is the classical Prandtl 
system. 

2.5 Quasi-hydrodynamic system for a viscous incompressible fluid 

In many particular cases of hydrodynamic flows we may neglect the density variation. Supposing 
that p and T are constant, from the equations 1)28(1 . (|29|) we obtain the system 



divu = divw, 



(32) 



du 1 ~ 1 

— + div(u <g> u) H — = _F -| (HvIIns + div \(w ® u) + (u ® w)l , (33) 

which is closed with respect to the unknown functions - the velocity u = u(x, t) and the pressure 
p = p(x, t). Here the vector 
w is defined by formula 

w = t{[u- V)u+ -Wp- Fj. 

We shall compute the tensor Iljvs using the expression 

Rns = ?][(V ® u) + (V ® w) T ]. 

The coefficient of dynamical viscosity r/ and the characteristic time r are considered to be given 
positive constants. Taking the formal limits in l|32l) ~ (|33() as r — > 0, we get the classical Navier- 
Stokes equations that describe the viscous non-compressible flows. The system (|32|I - H33|) is dis- 
sipative and possesses several explicit and physically reasonable solutions [3 - E]- The particular 
case of plane or spatial axially symmetric isothermal flows without external forces F = gives 



du x 1 d(y k u y ) 
dx y k dy 



dw x 1 d(y k w y ) 



dx 



y K 



On 



(34) 



du x 
~~dt 



d{u 2 x ) 1 d{y k u y u x ) 



+2 



dx 
d / du x 
dx V dx 
d(u x w x ) ^ 
dx 



V k dy 
]_d_ 

V k dy 
I d(y k u v w x ) 



1 dp 
p dx 
( du x du y 
dx 

1 d(y k u x w y ) 



y ^ dy 



dy 



dy 



(35) 



dUy 

~dT 



d{u x u y ) I d(y k u 2 y ) 1 dp 



d_ 

dx 
d(u x w v ) 



dx y k 
h ,( d u x du y 
\ dy dx 



dx 



d(u y w x ) 
dx 



dy 



pdy 



z a / k ou y 
y k dy \ dy 
2 d(y k u y Wy) 



dy 



2fcA 

y 



(36) 



where 



-{u x 



du x 
dx 



du x 



1 dp 
p dx 



w u 



'{u x 



dUy 

dx 



du y 
'~dy~ 



1 dp^ 
pdy i 



Here v = r// p is the coefficient of kinematic viscosity, the parameter k equals to zero in the plane 
case and equals to one in the axially symmetric one. The unknown variables are the components of 
the velocity u v = u y (x, y, t), u x = u x (x, y, t) with respect to the ortonornal local basis (e x , e y ) and 
the pressure p = p(x,y,t). The pressure field is defined using the already found fields of velocity 
and temperature by solving the Poisson equation: 



d 2 p Id , 
dx 2 y k dy 1 



d_ 

dx 



du x 
dx 



u 



dy 
V dy 



1 



du x 
dx 



1 d{y k u y ) 



y k dy 



yk Qy 



II 



f dUy dUy 



(37) 



This equation is the equivalent representation of (|34|l when r = const. 



3 Problem statement and computational algorithm 



Let us consider a plane two-dimensional incompressible flow in the channel of height H and of 
length L with small Mach numbers. The channel has a narrowing at the entrance section. The 
size of the narrowing is determined by the height of the step h. The scheme of the computational 
domain and the forming flow are demonstrated in Fig^ 
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Figure 1: Scheme of the computational domain 



We use the QHD system l|34[l - (|36[) with k — as the mathematical model. We transform this 
system into the dimensionless form, applying the relations 

x = xH, y = yH, u x =u x U , u y = u v U , p=ppUo, t=(tH)/U , Re=(U H)/u, 

where 

Uo = irrhj h My)dy 

is the flow velocity in the channel, averaged over the section, and uo(y) is the given velocity profile 
at the entrance section. We equip this dimensionless system with boundary conditions 
• the solid lower wall 



y = 0, < x < L/H, u x = u y = 0, — — = 0; 



y = 1, < x < L/H, u x = u y = 0, — — = 0; 



dp 
dy 

• the solid upper wall 

dp 
dy 

• the solid left wall 

dp 

x = 0, < y < h/H, u x = u v = 0, — = 0; 

ox 

• the inflow region at the left boundary 

dp 

x = 0, h/H<y<l, u x =u (y), u y = 0, — = const; 

• the right boundary 

du T du v 

x = LH, 0<y<l, —^ = —V=0, p = 0. 

dx dx 

Pressure boundary condition at the solid walls follows from the non-flow conditions for the 
velocity components and from the impermeability condition for the mass flow j m (|21|l . The pressure 
gradient at the channel entrance may be taken arbitrary. For example, it is possible to compute 
its values in the following way: we set the velocity profile at the channel entrance as the Poiscuillc 
parabola ^21) H3| : 

u (y) = ^^(l-y)(h/H-y). (38) 



The mass flow rate at the entrance section is computed according to the following formula 



/ / <><(0,y)-w x (0,y)]dy = -^(l-h/H) 3 ^-T(l-h/H)^. (39) 



h/H 



From P§)| we find 



dp 12J 



1 UT 



dx Re{l-h/HfV ' Re(l-h/H) 2 (J °' 

We chose the initial condition: u x = u y = 0. The pressure gradient at the initial moment was 
supposed to be constant all over the flow field. 

The dimensionless smoothing parameter r for laminar flows (|27|l was taken equal to 

7 Ma . , , Un _ c s H , , , 

r = + r , where Ma = — , i?e s = — (41) 

be Re s c s v 

- are the Mach number and the Reynolds number, derived from the speed of sound. For example, 
the air at normal temperature yields c s = 3.4 • 10 4 /, v = 0.15 2 /c, H = 10, Re s — 2 ■ 10 6 . For 
laminar flows we have Ma << 1. So for real flows the smoothing parameter proves to be small. 
We added to it the value To in order to compensate the difference scheme's antidiffusion and to 
provide stable computing. The value of to was chosen proportional to 1/ Re. 

The QHD equations are solved numerically using the algorithm, similar to the one described 
in ! - the explicit finite- difference scheme with second order of accuracy with respect to all 

spatial variables. Velocity and pressure values are defined in the same grid points. At each time 
step, the pressure field is calculated by using the velocity field, as a solution of Poisson equation 
(|37|) . which is also approximated with the second order space accuracy. The Poisson equation is 
solved by the preconditioned generalized conjugate gradient method. 

To present the numerical results, let us also introduce the stream-function, which is related to 
the solenoidal field u — w. These relations ^2] l°°k as follows: 

U »- W » = -^' U *- W * = ^dy- (42) 

The boundary conditions for the stream-function arc defined as follows - At the lower boundary of 
the computational domain and at the left wall we use the normalization -0 = 0, because there we 
have the impermeable boundary conditions. At the upper boundary the stream-function equals to 
the mass flow rate of the liquid. 



4 Numerical modeling of laminar flows 

For proper verification of the numerical method for back-step flow the problem described above 
has been solved with i?e=100, 200, 300, 400; h/H — 1/2. (From here on the Reynolds number is 
evaluated using the height of the step). The velocity profile at the entrance section represented 
the Poiseuille's parabola l|38|l . The dimensionless liquid mass flow rate J was taken equal to 1; it 
corresponded to the choice of the entrance pressure gradient in form of 



dp _ 96 
dx Re 



48r 
+ ~Re~ 



For the small values of t and the big values of Re we may suppose that 

dp 96 
dx Re 

The computed length of the separation zone behind the step was compared with data from [Hj. It 
was also defined from graphs presented in [Tnj . 

In [5] the Reynolds number was derived from the average flow velocity and from the height of 
the step. The entrance profile was also set in form of the Poiseuille's parabola. The mass flow 



rate J was taken equal to 1. The results, in particular, contain the length of the separation zone 
for H = 2h, Re(h)=l00, 200, 300. In ^U] the Reynolds number was derived from the value of 2h 
and the average entrance velocity. Graphic data concerning the length of the separation zone for 
50 < Re(2h) < 800 are presented here. 

The results, obtained by authors, are systematized in the tabled Here L is the dimensionless 
length of the computational domain, N x , N y are the numbers of mesh points in both directions, 
L s is the length of the separation zone, Nn er is the number of time steps till the conversion is 
achieved. The spatial mesh is uniform in both directions with equal widths h x = h y = 0.025. It is 
well known that the usage of equal widths h x and h y improves the accuracy of description of the 
separating flow. 

We have Re s ~ 10 6 in the described flows, so the value r = To in l|41|) was taken equal to 
To = 0.5/ Re. The time step St was equal to 10~ 4 for all variants of computation. 



Re(h) 


100 


200 


300 


400 


L 


7.5 


5.0 


7.5 


10 


N x X Ny 


300 x 40 


200 x 40 


300 x 40 


400 x 40 


T 


0.005 


0.0025 


0.00166 


0.00125 


Niter 


19800 


- 20000 


- 60000 


- 110000 


L s /h present comp 


5.0 


8.2 


10.1 


14.8 


L s /h |S], comput 


4.43 


7.5 


10.0 




L s /h 10 exp 


5.0 


8.5 


11.3 


14.2 


L s /h 1101 comput 


5.0 


8.3 


8.4 


7.8 



Table 1: Computations of laminar flows 
Computation stops when the condition 5p < 10~ 3 is satisfied. 

p n+l _ p n 

op = max , 

ot 

n is the time step number. 

In all variants the flow reaches the stationary regime. The length of the separation zone L s 
was defined by the location of the zero stream-function line. It is indicated with the accuracy 
0.2. Comparison of the results mentioned above with corresponding data from the Navier-Stokes 
simulation and with the experiments ^Oj demonstrates good agreement both in the length of the 
separation zone and in the picture of the flow in general. Mention, the good agreement for QHD 
and experimental results for Re = 400. An almost linear increase of the values of L s is observed 
in computations with the increasing number Re. 

For Re = 100 and 200 the process of flow relaxation consists of the appearence and further 
growth of a single vortex behind the step. For Re = 300 and 400 this process proves to be 
oscillatory and is accompanied with arising and separation of vortex-like formations, but, unlike 
the regimes with greater Reynolds numbers (they are considered in the second part of this paper), 
this oscillations fade and finally form a single stationary vortex behind the step. The isolines of 
the flow function tp, constructed according to 103), are demonstrated in Figs. El El They illustrate 
the process of flow relaxation in time for Re = 100 and 400. The isolines are placed equidistantly. 

With further increasing of the Reynolds number the stationary solution becomes unstable. 

The influence of the regularization parameter r and the conversion of the numerical solution 
was investigated for the variant with Re = 100. The value of t was additionally chosen equal to 
5 • 10~ 4 and 5 • 10~ 2 ; the time step St was changed proportionally. 

Besides the mesh described in the table, we also used another one - with twice as many nodes 
in both directions. The decreasing of the spatial mesh size by factor of two caused the analogous 
decreasing of the time step. It was shown that the length of the separation zone and the general 
picture of the flow practically doesn't depend neither of the value of regularization parameter r 
nor of the spatial mesh widths h x and h y . The increasing of r causes smoothing of the flow picture 



Figure 2: Stream functions for Re=100 




Figure 3: Stream functions for Re=400 

and allows us to increase the time step. Spatial mesh refinement gives a more detailed picture of 
the flow. 

We have studied the dependence of the solution on the pressure gradient at the entrance section 
with the average velocity and mass flow rate remaining constant. It was found out that the pressure 
gradient variation in the range from —96/ Re to —12/ 'Re practically doesn't influence the structure 
of the flow: at the distance around 0.5ft. from the entrance boundary the pressure adjusts to the 
existing liquid mass flow rate and practically doesn't depend on the initial gradient. 



5 Conclusion 



The present paper contains the phenomenological derivation of quasi-hydrodynamic equations. 
Two-dimensional mathematical model describing the viscous incompressible flow behind the back 
step is formulated and solved numerically. 

The computer simulation shows that the flows with small Reynolds numbers that correspond 
to the laminar regime, are stationary. The obtained regimes are in good agreement with the 
corresponding solutions of the Navier-Stokes system and with experimental data mentioned in 
literature. Oscillations that appear in the solutions describing the relaxation of laminar flows for 
moderate Reynolds numbers, fade with time. The final flow doesn't depend on the choice of the 
smoothing parameter r, which plays the role of regularizator in these computations. 

These results are in consistence with theoretical estimates f2L According to them, additional 
QHD-terms should be small in case of stationary flows and the solution of the QHD system is 
expected to be close to the solution of the Navier-Stokes system. Additional terms act as the 
regularizators and allow us to apply a relatively simple, stable and accurate numerical algorithm. 

The authors acknowledge Laboratoire D'Aerothcrmique du CNRS, Orleans, and personally Dr. 
J.-C. Lengrand, for the permanent support of this research. 
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